fitting.yb <- function(df, Col_x, Col_y, Col_group = NULL) {
  library(tidyverse)
  library(caret)
  library(ggplot2)
  
  if (!is.null(Col_group)) {
    df <- unite(df, "group", all_of(Col_group), sep = "_", remove = F)
    df$group <- as.factor(df$group)
    group <- unique(df$group)
  } else {
    df$group <- as.factor(1)
    group <- unique(df$group)
  }
  
  model <- data.frame(method =  rep("lm",6),
                      formula = c("y ~ x", "y ~ poly(x, 2, raw = TRUE)", "y ~ poly(x, 3, raw = TRUE)", 
                                  "y ~ poly(x, 4, raw = TRUE)", "y ~ poly(x, 5, raw = TRUE)", 
                                  "y ~ splines::bs(x, df = 3)"))
  result <- list()
  
  for (i in c(1:length(group))) {
    result[[i]] <- list()
    for (j in c(1:length(Col_y))) {
      df1 <- data.frame(x = df[df$group == group[i],Col_x],
                        y = df[df$group == group[i],Col_y[j]])
      df1$x <- as.numeric(as.character(df1$x))
      Model_performance <- data.frame()
      for (k in 1:6) {
        aa <- eval(parse(text = paste(model$method[k], "(", model$formula[k], ", data = df1)", sep = ""))) %>%
          summary()
        bb <- data.frame(aa$r.squared, pf(aa$fstatistic[1],aa$fstatistic[2],aa$fstatistic[3],lower.tail=F))
        Model_performance <- rbind(Model_performance, bb)
      }
      names(Model_performance) <- c("R2", "p")
      Model_performance$sig[Model_performance$p < 0.05] <- "*"
      row.names(Model_performance) <- c("line", "Polynomials^2", "Polynomials^3", "Polynomials^4",
                                        "Polynomials^5", "B-Spline")
      result[[i]][[j]] <- Model_performance
    }
  }
  return(result)
}